!=======================Resolve o sistema Ax=b=========================

PROGRAM exer05

!======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis========================

      INTEGER :: N, I, J
      REAL(4), ALLOCATABLE, DIMENSION(:,:) :: A, AUX, ADJ, MAT
      REAL(4), ALLOCATABLE, DIMENSION(:) :: b, x
      REAL(4) :: det

!==========Leitura de dados, abertura de arquivos e alocação===========

      WRITE(*,*) 'Entre com a dimensão da matriz A'
      READ(*,*) N

      OPEN(UNIT=10, FILE='equacao_in.dat', STATUS='old')
      OPEN(UNIT=11, FILE='equacao_out.dat', STATUS='unknown')

      ALLOCATE(A(N,N), ADJ(N,N), b(N), x(N))
      ALLOCATE(MAT(N,N), AUX(N-1,N-1))

      !Lendo do arquivo para a matriz/vetor
      READ(10,*) A
      READ(10,*) b

!==================Calcula a adjunta antes de transpor=================

      !Corre as linhas para gerar o valor
      DO I=1,N
            !Corre as colunas para gerar o valor
            DO J=1,N
                  !Calcula a matriz dos cofatores
                  CALL make_mat(A, I, J, N, AUX)
                  !Tira o determinante da matriz dos cofatores
                  CALL determinant(AUX, N-1, det)
                  !Calcula o elemento I,J de MAT
                  MAT(I,J) = ((-1)**(I+J)) * det
            END DO
      END DO

!======================Obtém a matriz inversa de A=====================

      !Calcula o determinante para depois dividir pela adjunta
      CALL determinant(A, N, det)
      !Calcula a adjunta (que é a transposta de MAT)
      ADJ = TRANSPOSE(MAT)
      !Calcula a inversa que é a adjunta dividia pelo determinante
      A = (1/det)*ADJ

!=====Obtém x multiplicando a inversa por b e escreve o resultado======

      !Cálculo de x (x = b A^-1)
      x = MATMUL(b,A)
      !Escreve o resultado no arquivo de saída
      WRITE(11,*) x

!======================================================================

END PROGRAM exer05

!======================================================================


!==============================Subrotinas==============================

!========================Faz a matriz dos cofatores====================
!Recebe a matriz de entrada, a linha e a coluna a serem ignoradas,
!a dimensão N da matriz, e o vetor AUX que é o retorno da função

SUBROUTINE make_mat (MAT_IN, I, J, N, AUX)

      IMPLICIT NONE

      INTEGER :: I, J, N, K, L, linha, coluna
      REAL(4), DIMENSION(N,N) :: MAT_IN
      REAL(4), DIMENSION(N-1,N-1) :: AUX

      !Inicialica as variáveis que dizem onde entra o elemento na
      !nova matriz.
      linha = 1
      coluna = 1

      !Loop que percorre as linhas da MAT_IN
      DO K=1,N
            !Elimina a linha em que estamos
            IF (K /= I) THEN
                  coluna = 1
                  !Loop que percorre as colunas da MAT_IN
                  DO L=1,N
                        !Elimina a coluna em que estamos
                        IF (L /= J) THEN
                              !Adiciona elemento a nova matriz
                              AUX(linha,coluna) = MAT_IN(K,L)
                              !Incrementa a coluna para saber onde
                              !escrever depois
                              coluna = coluna + 1
                        END IF
                  END DO
                  !Incrementa a linha para saber onde escrever depois
                  linha = linha + 1
            END IF
      END DO

END SUBROUTINE make_mat



!=======Subrotina que calcula o determinante de uma matriz NxN=========

SUBROUTINE determinant(M_IN,N,det)

      IMPLICIT NONE

      INTEGER :: N, fat, I, J
      INTEGER, DIMENSION(N+1) :: vet
      REAL(4), DIMENSION(N,N) :: M_IN
      REAL(4) :: det, termo

      !Abre arquivo inicial para gerar as permutações
      OPEN(UNIT=70, FILE='perm_in.dat', STATUS='unknown')
      !Escreve '1 1' no arquivo para poder gerar as permutações
      WRITE(70,*) '1 1'
      !Fecha arquivo gerado
      CLOSE(UNIT=70)

      !Chama as permutações (necessário para o cálculo)
      CALL make_perm_out(N)

      !Abre o arquivo onde se encontram as permutações
      OPEN(UNIT=40, FILE='perm_out.dat', STATUS='old')

      !Inicializa det
      det = 0.0

      !Calculo do fatorial que é usado como limite no proximo loop
      fat = 1
      DO I=1,N
            fat = fat*I
      END DO

      !Vai lendo linha por linha do arquivo as permutações P_i e sua
      !paridade para com a ajuda delas, calcular o determinante por
      !uma série
      DO I=1,fat
            termo = 1.0
            !Lê a permutação P_i do arquivo
            READ(40,*) vet
            DO J=1,N
                  !Calcula o termo a ser somado na série
                  termo = termo*M_IN(J,vet(J))
            END DO
            !Soma o termo na série (sinal depende da paridade)
            det = det + vet(N+1)*termo
      END DO

      CLOSE(UNIT=40)
      CALL SYSTEM('rm -f perm_out.dat')

END SUBROUTINE determinant



!=====Subrotina que gera as permutações para poder calcular o det======

SUBROUTINE make_perm_out (N)

      IMPLICIT NONE

      INTEGER :: N, I, J, fat
      INTEGER, ALLOCATABLE, DIMENSION(:) :: leitura

      DO  I=1,(N-1)
            OPEN(UNIT=30, FILE='perm_in.dat', STATUS='old')
            OPEN(UNIT=20, FILE='perm_out.dat', STATUS='new', POSITION='append')

            ALLOCATE(leitura(I+1))

            !Calcula o fatorial que será usado como limite no proximo loop
            fat = 1
            DO J=1,I
                  fat = fat*J
            END DO
            !Lê linha por linha do arquivo de entrada e permuta a partir de
            !cada linha. Chama a subrotina para cada linha
            DO J=1,fat
                  READ(30,*) leitura
                  CALL permut(leitura,I)
            END DO
            
            CLOSE(UNIT=30)
            CLOSE(UNIT=20)

            DEALLOCATE(leitura)

            CALL SYSTEM ('mv perm_out.dat perm_in.dat')
      END DO

      !Muda o nome do arquivo para ficar como saída
      CALL SYSTEM ('mv perm_in.dat perm_out.dat')

END SUBROUTINE make_perm_out



!=================Subrotina usada pelo make_perm_out===================
!Subrotina recebe o vetor no qual incluirá o próximo número e fará
!as permutações. Também recebe N (número inicial de valores fora a
!paridade)
SUBROUTINE permut (vet,N)

      IMPLICIT NONE

      INTEGER :: N, I, J, flag
      INTEGER, DIMENSION(N+1) :: vet
      INTEGER, DIMENSION(N+2) :: temp

      DO J=1,(N+1)
            !Inicializa o vetor temporário
            temp = 0
            !Verifica a paridade
            IF (vet(N+1) == 1) THEN
                  temp(N+2) = (-1)**(J+1)
            ELSE
                  temp(N+2) = (-1)**J
            END IF
            !Coloca o valor de N+1 na posição que depende de J
            temp(N+2-J) = (N+1)
            flag = 0

            !Loop que faz as permutações
            DO I=1,N
                  IF (temp(I) /= (N+1) .AND. flag == 0) THEN
                        !Para os valores diferentes e anteriores a N+1
                        temp(I) = vet(I)
                  ELSE
                        !Para os valores após o N+1
                        temp(I+1) = vet(I)
                        flag = 1
                  END IF
            END DO

            !Escreve no arquivo de saída as permutações geradas nesse loop
            WRITE(20,*) temp
      END DO

END SUBROUTINE permut



!======================================================================
